/***************************************************************************
 *   Copyright (C) 2006 by Harish Kukreja                                  *
 *   harish.kukreja@gmail.com                                              *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

//
// C++ Interface: physiksvector
//
// Description: Wrapper class around GSL's gsl_vector. Keeps SWIG happy and developers sane.
//
//
// Author: Harish Kukreja <harish.kukreja@gmail.com>, (C) 2006
//
// Copyright: See COPYING file that comes with this distribution
//
//

#ifndef _PHYSIKSVECTOR_H_
#define _PHYSIKSVECTOR_H_

#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>

namespace PhysiksBase {

/** \class Vector
 * A simple wrapper class around gsl/gsl_vector. Without this, one would have an unpleasant
 * blend of C-style code alongside C++ code. Note that it is work in progress, with new
 * functions added when they are needed. */
class Vector {
    public:
        /** \brief  Constructor 1.
         * \param size The dimension of the vector
         * \param init the initialisation value of all vector elements */
        Vector(unsigned int size, double init = 0.0);

        /** \brief  Constructor 2.
         * \param init the initialisation gsl_vector */
        Vector(const gsl_vector *init);

        /** \brief  Copy Constructor.
         * \param init the reference Vector */
        Vector(const Vector &init);

        virtual ~Vector();


        /** \brief  In place scalar multiplication (i.e. with side effect).
         * \param fac the scaling factor */
        void scale(double fac);

        /** \brief  Set the vector to the zero vector. */
        void set_zero();

        /** \brief  Set the value of a coordinate. No bounds checking.
         * \param index the index of the coordinate (0-based)
         * \param value the value to change the coordinate to */
        void set(unsigned int index, double value);

        /** \brief Sets the coordinates of this vector to those of the parameter.
         * \param copy the template vector, of which coordinates are to be copied over  */
        void set(const Vector &copy);
        
        /** \brief  The Euclidean norm (magnitude) of the Vector.
         * \return the euclidean norm. */
        double norm() const;

        /** \brief  A unit vector in the same direction as the calling vector.
         * \return a unit vector in the same direction as the calling vector */
        Vector normalize() const;
        
        /** \brief Effectively "clones" the current vector.
         * \return  an identical copy of this vector    */
        Vector duplicate() const;

        /** \brief  Returns the dimension of the vector.
         * \return calling vector's dimension */
        unsigned int dim() const;

        /** \brief  Get accessor for an arbitrary coordinate. Bounds checking offered.
         * \param index The coordinate's index (0-based).
         * \return the value of the respective coordinate. 0 if index is out of bounds */
        double get(unsigned int index) const;
        
        /** \brief  Vector cross product. Does not check dimensions of cross-product operands, so
         *          vectors of dimension other than 3 might yield non-sensical results.
         * \param   rhs the 2nd operand of the cross product.
         * \return  the resultant vector, orthogonal to both the caller and rhs. */
        Vector cross(const Vector &rhs) const;


        /** Vector addition overload. */
        friend Vector operator+(const Vector &lhs, const Vector &rhs);
        /** Vector addition & assignment overload. */
        friend void operator+=(Vector &lhs, const Vector &rhs);
        /** Vector subtraction & assignment overload. */
        friend void operator-=(Vector &lhs, const Vector &rhs);
        /** Vector subtraction overload. */
        friend Vector operator-(const Vector &lhs, const Vector &rhs);
        /** Vector negation overload. */
        friend Vector operator-(const Vector &rhs);
        /** Vector standard inner product (dot product) overload. */
        friend double operator*(const Vector &lhs, const Vector &rhs);
        /** Vector scalar multiplication overload. */
        friend Vector operator*(const Vector &lhs, double fac);
        /** Vector scalar multiplication overload. */
        friend Vector operator*(double fac, const Vector &rhs);
        /** Vector indexing (random coordinate access) function. No bounds checking. */
        double operator[](unsigned int k) const;

        /** \brief  Returns a pointer to the gsl_vector core of the calling Vector.
         * PLEASE DO NOT USE THIS. It exists solely for (rare) internal usage. */
        gsl_vector *core() const;

    private:
        /** The dimension of the Vector. This should be immutable, no? */
        unsigned int dimension;
        
        /** The gsl_vector core. */
        gsl_vector *cont;
};

}   // End namespace PhysiksBase.

#endif
